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ABSTRACT 

We investigate whether any multi-planet systems among i^^ep/er candidates (2011 February release) 
can harbor additional terrestrial-mass planets or smaller bodies. We apply the "packed planetary 
systems" hypothesis that suggests all planetary systems are filled to capacity, and use a Hill stability 
criterion to identify eight 2-planet systems with significant gaps between the innermost and outermost 
planets. For each of these systems, we perform long-term numerical integrations of 10^ years to 
investigate the stability of 4000—8000 test particles injected into the gaps. We map out stability regions 
in orbital parameter space, and therefore quantify the ranges of semi-major axes and eccentricities of 
stable particles. Strong mean-motion resonances can add additional regions of stability in otherwise 
unstable parameter space. We derive simple expressions for the extent of the stability regions, which 
is related to quantities such as the dynamical spacing A, the separation between two planets in units 
of their mutual Hill radii. Our results suggest that planets with separation A < 10 are unlikely to host 
extensive stability regions, and that about 95 out of a total of 115 two-planet systems in the Kepler 
sample may have sizeable stability regions. We predict that Kepler candidate systems including KOI 
433, KOI 72/Kepler-lO, KOI 555, KOI 1596, KOI 904, KOI 223, KOI 1590, and KOI 139 can harbor 
additional planets or low-mass bodies between the inner and outer detected planets. These predicted 
planets may be detected by future observations. 

Subject headings: planetary systems - methods: numerical - planets and satellites: dynamical evo- 
lution and stability - stars: individual (KOI 433, KOI 72/Kepler-lO, KOI 555, 
KOI 1596, KOI 904, KOI 223, KOI 1590, KOI 139) 



1. INTRODUCTION 



Early studies of extrasolar planetary systems showed 
residual velocity tr ends in Keplerian orbit fits to radial 



veloc i ty data (e.g., |Marcy fc Butler||1 998; But ler et al. 
19981 IMarcy et al ll999| |Vogt et al. 2000( [Fisch er et al. 
200T|, suggesting that the se s yste ms may hos t addi- 
tional, undetected planets. [Fischer et al.| ( |200l| ) noted 
that about half of the stars in their sample of 12 systems 
showed residual trends greater than the expected scat- 
ter due to measurement uncertainties and stellar noise. 
Most of these systems were later confirmed to harbor 
additional planet (s). 

In more recent years, the study and prediction of undis- 
covered planets have been aided by long-term N-body 
simulations. These numerical investigations searched for 
stability zones in multi-planet systems by integrating 
hundreds to thousands of test bodies, which were injected 
into empty regions be tween known plane ts (e.g., Menou 
fc Ta bachnik 2003; B arnes fc Raymond|[2Q 04 ; Raymond 
F^rnes..2005, .Ji et al.||2005[ [Rivera fc 
20071 IRaymond et al.||2008D. F( 



74 



iver a &: Hagh ighipour 
or example, a putative 



Saturn-mass planet in HP t4156 was first predicted by 
[Raymond fc Barnes| (|2005[) through numerical simula- 
tions that showed a stable region between planet s b and 

' ' ' ' ' ' ~ ' ( [2008) , 



The planet was later discovered by Bean et al. 



although there have been questions about the validity of 



the d etection (Wittenmyer et al. 2009 Meschiari et al. 



2011[ ) . This prediction was motivated by the ''packed 
planetary systems" (PPS) hypothesis. 
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The PPS hypothesis is the idea that planetary sys- 
tems are formed "dynamically full" and filled to capac- 
ity, and any additi onal planets will cause the systems to 
be unstable (e. g., Barnes fc Ra ymond 2004; Raymond 

Barnesll2005l [Raymond et al.||2006^ ,Barnes fc Green- 



berg 2007|. Consequently, planetary systems with sta- 



ble stability zones between the innermost and outermost 
planets are suggestive of additional, yet-undetected plan- 
ets. Reasons for the non-detections of hypothetical plan- 
ets include lack of sufficient data, such as non-transiting 
planets that require more data to detect them via tran- 
sit timing variations, or planetary masses/radii that are 
below detection limits. The orbital properties of pre- 
dicted planets can be identified through long-term nu- 
merical simulations. Support for the PPS hypothesis 
comes from early observations of packed multi-planet 
systems that led to this hypothesis (e.g ., iButler et al. 



19991 

et al 



Marcy et alT]|2001a|b| [Fischer et al. 2002; [Mayor 



2004), apparent consistency between t he planet 



planet scattering model and packed systems ([Raymond 



et al. 2009 L the re markably dense and packed Kepler- 



11 system~i[Lissauer et al. 2011a), theoretical work (e.g.^ 
Chambers et al.|l§5 6 ; Smith fc Lissauer 2009 J, and other 
investigations (e.g., .Rivera fc Lissauer|2000f [Gozdziewski 



Migaszewski 2006f~ 
In the present study, we apply the PPS hypothesis to 
multi-planet candidate systems discovered by the Kepler 
team du ring the mission's fir st four and a half months 



of data (Borucki et al. 



2011). The Kepler mission is a 



transit survey designed to search for Earth-sized plan- 



ets (Boruc ki et al.[[2Q10j [Koch et al.||2QTQ| [Jenkins et al 
2010; Cal dwell et al.|2010l ), and is sensitive to terrestrial 



class and larger planets located at a large range of sep 



arations from their host star. Kepler can detect mul- 
tiple transiting systems for densely-packed planets with 
nearly coplanar configurations or with serendipitous ge- 
ometric alignment, and the dynamics and statistics of 
Kepler multi-planet systems are providing a wealth of 
inforrnation about planetary systems (e.g., Steffen et al.„ 
2Q1Q[ [Latham et al.||2QlT| [Lissauer et al.|2QTTa|F| r 



Uiven that planetary systems have been discovered 
with densely packed planets, we seek to test the PPS 
hypothesis and predict additional planets in Kepler can- 
didate multi-planet systems. In Section [2j we discuss 
Kepler'^s sample of multi-planet systems as well as our 
methodology for evaluating each planetary system's level 
of dynamical packing. We also explain our methods for 
running numerical simulations and our choice of initial 
conditions. In Section |3) we present the results from 
long-term numerical integrations and illustrate them us- 
ing stability maps. Section [4] discusses the dynamical 
interpretation of our work, in particular the relationship 
between planetary spacing and the extent of an inter- 
planet stability region. We then summarize the restric- 
tions and scope of our study (Section [5| and state our 
conclusions (Section [6|. 

2. DATA AND METHODS 

Based on publicly- available Kepler data covering the 
first four and a half months of observations, about one- 
third of '^1200 transiting planet candidat es a r e hosted 
in multi-plan et systems (Borucki et al. 2 Q11[ [Lissauer 



et al. 



20TTb|. 



These multi-planet systems mclude 115 
2 transiting planets, 45 systems with 3 
8 systems with 4 transiting planets. 



systems with 
transiting planets 
1 system with 5 transiting planets, and 1 system with 6 
transiting planets. Most Kepler candidate planets have 
not been validated and are therefore "Kepler Objects 
of Interest" (KOI) and assigned a number. Candidates 
in multi-planet systems have a smaller probability than 
single-pl anet candidates of being an astrophysical false 
positive d Ragozzine fc Holman |2QlQ Latham et al.'2011^ 
[Lissauer et al. 2011a, 2012). Moreover, all of these candi- 
date multi-planet systems are stable over long-term inte- 
grations (Lissauer et al. 20 11b). For the remainder of this 
paper, we refer to all candidate planets and systems by 
dropping the adjective "candidate. " All of these Kepler 
multi-planet systems presented by [Borucki et al. (2011) 
are examined using the analytical method described be- 
low in Section \2A\ 

To discern the extent of packing in Kepler multi-planet 
systems, we defi ne two types of stability as outlined by 



Gladman (1993). Fulfillment or over- fulfillment of these 
stability criteria, meaning that the considered planetary 
system is not on the verge of instability, can imply the 
presence of additional planets according to the PPS hy- 
pothesis. First, Hill stability requires that a system's 
ordering of planets (in terms of distance from the star) 
remains constant. This means that close approaches are 
forbidden and planet-crossing is not allowed for all time, 
but the outer planet may be unbound and still be Hill 
stable. The second type of stability is Lagrange stability, 
which is a stricter definition than Hill stability. Lagrange 
stability requires not only the conservation of the order- 
ing of planets, but also that they remain bound to the 
star for all time. Hill stability can be mathematically 
examined for two-planet, non-resonant systems, and La- 



grange stability is typically examined through numerical 
simulations. Hill stability can be a reasonable p redictor 
of Lagrange stability ( |Barnes fc Gree nbergl2 006[ ) . In the 
next two subsections, we exa mine Hill stability through 
analytical methods (Section [2.1[ ) and Lagrange stabil- 
ity through N-body integrations (Section 2.2) for Kepler 
multi-planet systems. 

2.1. Analytical Method 

We calculate the Hill stability criterion of each adjacent 
planet pair in ii^e^/er multi-planet systems. W e follow the 
notation given in Barnes fc Greenberg (2006) by referring 
to the relevant quantities as /j an3"^rit* 



-2(M*+Mi+M2) 



=1 + 3^/^ 



M1M2 



M*/^(Mi 
MiM2(llMi -i 



+ M2)4/3 
7M2) 



(1) 



(2) 



3M*(Mi + M2)2 

where is the mass of the star. Mi and M2 are the 
masses of the planets where Mi > M2, L and E are 
the total orbital angular momentum and energ y of the 
system, and G is the gravitati onal constant ( Marchal 
fc Bozisj|l982l [Gladman||l993| [Veras fc Armita^jpnUll r 
Iwo-planet, non-resonant systems with /i^//icrit ^ 1 are 
considered Hill stable. For systems with additional plan- 
ets and/or in resonance, stability needs to be investigated 
numerically. A system that does not fulfill the Hill sta- 
bility criterion has unknown Hill stabilityjit may or may 
not be Hill stable. In Equations ([T| and ([2|, /3crit is only 
a function of masses and (3 is a function of masses as 
well as semi-major axes and eccentricities; evidently, for 
a given set of masses, there are stability boundaries in 
orbital parameter space. 

We calculate /3//3crit for each adjacent planet pair in 
Ke pler multi-planet syst ems released as of February 2011 
by Borucki et al.[ ( [2011 ) . We find that almost all of the 
adjacent planet pairs h ave /3//3crit values greater than 1. 
Raymond et al. (2009) discussed that planet pairs with 
values of fJ / /:/crit greater than ~1.5— 2 are probably capa- 
ble of harboring additional planet (s) with a semi-major 
axis in between those of the existing planets. We find 
eight iTep/er systems, all of which are two-planet systems, 
with j3 /j3cT\t values of 1.5 or greater (Table [T]). None of 
the 3-planet, 4-planet, 5-planet, and 6-planet systems 
have any adjacent planet pairs with /3//3crit values of at 
least 1.5. In the next section, we place test particles in 
each of these eight systems to determine their zones of 
stability. 

Another useful criterion for evaluating stability is the 
dynamical spacing A between two planets, i.e., the dif- 
ference between their semi- major axes expressed in units 
of their mutual Hill radius. 



A = 



a2 - ai 



where Rhi 2 is the mutual Hill radius defined as 



R 



Hl,2 



Ml +M2 
3M* 



1/3 



ai + a2 



(3) 



(4) 
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Table 1 

Identified Kepler Systems with /5//3crit > 1-5 



KOI 


M* (Mo) 


Ml (Me) 


Ri (i?©) 


ai (AU) 


Pi (days) 


Ms (Me) 


R2 (i^e) 


as (AU) 


P2 (days) 


/3//3crit 


A 


433 


1.01 


37.38 


5.80 


0.050 


4.030 


209.81 


13.40 


0.935 


328.240 


2.861 


28.7 


72t 


1.03 


1.72 


1.30 


0.018 


0.837 


5.56 


2.30 


0.252 


45.295 


2.781 


90.3 


555 


0.95 


2.31 


1.50 


0.046 


3.702 


5.56 


2.30 


0.376 


86.496 


2.031 


77.3 


1596 


0.87 


5.56 


2.30 


0.061 


5.924 


13.21 


3.50 


0.416 


105.355 


1.817 


53.4 


904 


0.69 


4.61 


2.10 


0.029 


2.211 


9.61 


3.00 


0.159 


27.939 


1.624 


50.4 


223 


0.92 


7.74 


2.70 


0.041 


3.177 


6.07 


2.40 


0.226 


41.008 


1.621 


56.2 


1590 


0.88 


3.75 


1.90 


0.033 


2.356 


8.34 


2.80 


0.163 


25.780 


1.527 


55.4 


139 


1.07 


1.46 


1.20 


0.045 


3.342 


36.07 


5.70 


0.741 


224.794 


1.508 


54.1 



Eight Kepler multi-planet systems are identified with (3/(3crit values greater than 1.5. We list their KOI ("Kepler Objects of 
Interest") identifier, stellar mass M*, planetary masses Mi and M2, planetary radii Ri and R2, semi-major axes ai and as, 
orbital periods Pi and P2, /3//3crit value (listed in descending order), and dynamical spacing criterion A. The inner planet has the 
subscript 1 an d the outer plan et has the subscript 2. Stellar mass and planetary parameters (size, semi-major axis, and period) 
are taken from|Borucki et al.'f2011), and we derived the planetary masses using a power law: Mi = {Ri / R^)'^-^^ where the 

subscript i represents planet 1 or 2 (Lissauer et al. 2011^. 

t KOI 72 is a confirmed system and is also known as Kepler- 10 ( [Batalha et al.|2011| ). 



(e. 

scripts 



Gladman 
1 and 



1993[ [Chambers et al |[T996l ). Here, sub 



spectively. 



denote the inner and outer planets, re 



2.2. Numerical Method 



The previous section identified eight Kepler systems 
with /3//3crit > 1-5 (Table [T), which suggests that these 
systems are most likely tonave gaps between adjacent 
planets that may contain additional planet (s). For these 
eight systems, we numerically explore their regions of 
Lagrange stability to determine zones in orbital element 
space that can harbor additional, undetected planets 
that are stable. We use a hybrid symplectic/Bulirsch- 
Stoer alg orithm from an N-body integration package. 
Mercury ( Chambers|1999 ) , with a timestep that sampled 
1/20*^ of the innermost planet's orbit. 

For each of the eight identified Kepler systems, our 
simulations include the star and its two detected planets 
as well as 4000—8000 massles^ test particles placed in 
between the locations of the inner and outer planets. We 
do not assign a common number of test particles to each 
system for computational cost reasons. The overall cost 
of the integration is a function of the timestep and of the 
number of test particles. Systems that have inner plan- 
ets with shorter orbital periods (and therefore shorter 
timesteps) are assigned fewer test particles so that the 
integrations may finish within a reasonable amount of 
time. In total, we integrate '^8000 test particles per sys- 
tem except for KOI 72 (--4000 test particles) and KOI 
904 (~6000 test particles). Each system is integrated for 
lO*" years, and test particles that survive the length of 
the integration are considered stable test particles. 

For each of the eight identified Kepler systems, ini- 
tial conditions for the star and its two detected plan- 
ets are given in Table [l] These initial conditions in- 

^ It would be ideal to perform detailed integrations of numerous 
test bodies with nonzero masses, distributed with varying distances 
and velocities from the star to sample all possible orbits. However, 
this process would be very computationally costly. Since we are in- 
vestigating the stability of terrestrial-mass planets or smaller bod- 
ies, we approximate such objects as massless test particles in our 
simulations. These test particle approximations have been simi- 
larly adopted in previous studies (i.e., |Rivera Haghighipour|2007| 
and references therein). 



elude the star's mass and the known planets' masses, 
radii, semi-major axes, and orbital periods. The other 
orbital elements of the planets-eccentricity, inclination, 
argument of pericenter, longitude of the ascending node, 
and mean anomaly-are currently unknown; we assume 
circular and coplanar orbits and assign a random mean 
anomaly for the planets. The coplanar assumption is 
supported by the fact that these are all transiting plan- 
ets; the larger the mutual inclination between the plan- 
ets' orbital planes, the smaller the probability t hat they 
all transit the star (Ragozzine & Holman |2010|). Circu- 



lar and coplanar orbits have the least angular momentum 
deficit an d therefore a re most likely to be stable config- 
urations ( Laskar|[l997 ). 

Initial conditions tor test particles are as follows. In- 
clinations i are drawn from a uniform distribution ( 0° < 
i < 5°). Previous work by Lissauer et al. (2011b) ini- 
tially suggested that Kepler multi-planet systems have 
low relative inclinations with a mean of ^5°, but that 
number was revised to ^10° as our paper was undergo- 
ing revisions. Semi-major axis a and eccentricity e are 
initially drawn from a uniform distribution (ai < a < a2; 
< e < 1) for the first 1000 particles for each system. 
Subsequent integrations of additional test particles were 
randomly inserted into semi-major axis and eccentricity 
bins that had few or no particles, by filling up bins with 
lower eccentricity first. This ensured better coverage of 
semi-major axis and eccentricity space. All other orbital 
elements (argument of pericenter, longitude of the as- 
cending node, and mean anomaly) are drawn randomly 
from a uniform distribution between 0° and 360°. 

This procedure is performed for each of the eight Ke- 
pler systems identified with /3//3crit > 1-5 (shown in Table 
n]). For each system, we record each test particle's start- 
mg orbital elements and whether it became unstable or 
remained stable during the duration of the integration. 
Instability can be due to ejection or collision of the test 
particle with another body. 

3. RESULTS 

In this section, we describe the results stemming from 
long-term N-body integrations of eight Kepler systems 
(Table [1]). All of these multi-planet systems have two 
known planets, and were identified in Section [21!] as po- 
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Figure 1. Stability maps for KOI 433, KOI 72/Kepler-lO, KOI 555, and KOI 1596. Test particles, shown as circles, 
are displayed at their starting values of semi-major axis and eccentricity. Filled blue circles are test particles that 
survived the integration length of 10^ years, and unfilled gray circles are test particles that did not remain stable in 
that time. Black curves show the boundaries dividing planet crossing and non-planet crossing orbits. Vertical black 
lines represent the locations of first-order mean- mot ion resonances with the outer planet. The inner and outer planets 
are located at the left and right edges of the plot, respectively. 



tentially capable of containing an additional planet in 
the regions between the inner and outer planets. For 
each of these systems, we quantify their zones of sta- 
bility and instability in semi-major axis and eccentricity 
space by illustrating them in stability maps plotted in 
Figures [l]-[^ 

Our stability maps and results indicate that each of 
these planetary systems are capable of harboring a sta- 
ble low-mass body for up to 10^ years in the interme- 
diate zone between the known inner and outer planets. 
We discover broad stable regions in each planetary sys- 
tem, which appear as mountain-shaped regions in Fig- 
ures [l]-[2] We also find additional regions of stability 
outside the mountain regions, where test particles can 
have stable orbits due to mean-motion resonances with 
the inner and outer planets. Strong first-order resonances 
with the outer planet are marked in Figures [T]-[2j As for 
instabilities, the majority (typically ~90%) of unstable 
test particles were unstable within the first 10^ years. 



Stable test particles do not show much movement in 
semi-major axis and eccentricity over the course of an 
integration. As a result, the plots shown in Figures [l]-[2] 
only show the starting locations of test particles. We 
quantify the motion of test particles in orbital element 
space by computing the median of the absolute values of 
the differences between initial and final values of semi- 
major axis and eccentricity. The median semi-major axis 
differences range from ~ 3.8 x 10~^ AU to ~ 2.3 x 10~^ 
AU, and the median eccentricity differences range from 
- 2.4 X 10-"^ to - 8.5 X 10-"^. The largest differences 
in semi-major axis and eccentricity for stable particles 
are commonly due to particles placed near the edge of 
the stability region that became scattered off to another 
part of the stability region, or particles originally not 
in the stability region that became scattered to an or- 
bit with a final semi-major axis greater than the outer 
planet's semi- major axis and typically accompanied by 
an increase in eccentricity. 
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Figure 2. Same as Figure [TJ except we show stability plots for KOI 904, KOI 223, KOI 1590, and KOI 139. 



Mean motion resonances can act as additional reser- 
voirs of stability outside of the mountain-shaped region. 
Strong first-order resonances are plotted in Figures [l]-[2] 
to provide examples of stable test particles in resonances 
outside of the stability region. Many more first-order and 
higher-order resonances exist, forming a thicket of reso- 
nance locations that are not drawn to reduce confusion. 
We find that the majority of stable test particles outside 
the mountain are located in resonant or near-resonant 
periods with that of the inner or outer planet. Test par- 
ticles placed in planet-crossing orbits (above the black 
curves drawn in Figures [l]-[2| can be stable if placed in 
such resonances, which protect the test particles from 
close encounters with the planets. The role of mean mo- 
tion resonances in the stability of test particles and plan- 



ets has also been previously explored (e.g., see Rivera & 



H aghighipour||2007| |Barnes fc Gree nberg^'2007') . 

Our results show that massless test particles can stably 
orbit in these stability regions for up to 10^ years, and we 
suggest that these stability results can be extended from 
massless particles to Earth-mass planets. Spot checks 
performed for KOI 1596, which has a moderate (for this 
sample) /3//3crit of ~ 1.8 17, indicate that an Earth-mass 
planet with a semi-major axis in the middle of the sta- 



bility region and with an eccentricity of zero is stable for 
at least 10^ years. We have also tested scenarios where 
we increased the mass of the inserted planet up to a few 
Earth masses, as well as cases where we inserted two, 
three, and four evenly-spaced. Earth-mass planets with 
zero eccentricities in the main stability region. These in- 
tegrations all proved to be stable for up to 10^ years in 
our tests for KOI 1596. As a result, it is likely that the 
stability zones identified using massless test particles are 
applicable to Earth-mass bodies, and that these stability 
zones can potentially contain more than one Earth-mass 
planet. 

We briefly compare our numerical results with analyti- 
cal expectations. Our stability results based on this sam- 
ple of Kepler systems indicate that two-planet systems 
meeting the analytical threshold /3//3crit > 1.5 are consis- 
tent with the idea that they can hold additional planet (s) 
in intermediate separations from their host star. All eight 
systems investigated here had planet pairs with /3//3crit > 
1.5 and were numerically found to be "unpacked," which 
supports previous work suggesting that additional plan- 
ets are expected to be stable in systems with ^ / ^cvit > 
1.5—2 ([ Raymond et al.|2009 ) . Since all eight Kepler sys- 
tems we investigated with 13 / ^ exit > 1.5 are unpacked, we 
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expect that there may also be additional systems with 
P/Pcvit less than 1.5 that are also unpacked (see next 
section). 

4. PLANETARY SPACING DETERMINES EXTENT OF 
STABILITY REGION 

We now describe in greater detail the shapes and sizes 
of the mountain-shaped stability regions observed in Fig- 
ures [T]-[2] In particular, we discuss relationships between 
the spacing between two planets and the extent of the 
stability region in-between the planets. 

The stable regions in each planetary system include 
a mountain-shaped stability peninsula as well as nar- 
row strips of stability due to mean-motion resonances 
(e.g., see Figures IT]-[2| . The large mountain-shaped sta- 
bility region has a shape common to all of the plane- 
tary systems, because it is sculpted on the left and right 
flanks by specific semi-major axis and eccentricity values 
that delineate planet-crossing orbits. Mathematically, 
the mountain's left slope is shaped by ai — a(l — e), 
where a\ is the inner planet's semi-major axis. The 
mountain's right slope is shaped by a2 = a(l + e), where 
a2 represents the outer planet's semi-major axis. These 
orbit-crossing boundaries are shown as black curves in 
each stability plot in Figures [l]-[2] The actual stability 
boundaries (left and right fianks of the mountain- shaped 
stability region) do not extend all the way to the black 
curves. This is explained by close approach effects: test 
particles that are not initially on planet-crossing orbits 
can become unstable if they make sufficiently close ap- 
proaches to the existing planets. The critical distance 
from the planet-crossing boundary at which this can oc- 
cur is similar to the half- width of a planet's "feeding 
zone" in which planetesimals may impact the planet, 
which can be estimated at about '^2.3 Hill radii for circu- 
lar orbits (i.e., |Greenberg et al.|1991| ). The results of our 
simulations show similar distances between the planet- 
crossing curve and the actual slope of the mountain. 

The maximum height of each mountain-shaped sta- 
bility region is constrained by the semi-major axes of 
the inner and outer planets. The maximum eccentricity 
allowed is determined by the intersection of the orbit- 
crossing boundaries, a\ = a(l — e) and a2 = a{l + e), 
and serves as an upper limit to the maximum possible 
height of the stability mountain. Since we assume that 
the two known planets in each system have orbital ec- 
centricities of zero, the intersection of the curves occurs 
at a semi-major axis of (ai + a2)/2 and an eccentricity 
of 



= 1 



2ai 



(5) 



which is the maximum possible eccentricity emax of the 
stability mountain. As evident in Figures [l|-[2j the actual 
peak of the stability region is not the same as the Cmax- 
The actual height or peak of each mountain-shaped 
stability region can be computed as follows. Consider a 
test particle located between the inner and outer planets. 
In order to remain stable, this particle cannot enter a 
zone of dynamical inffuence surrounding each planet. We 
measure this exclusion zone as a certain number q of 
mn radii Rhz, where Rhz = (Mi/(3M*))^/^ai and the 
subscript z = 1,2 refers to the inner and outer planets, 
respectively. Therefore, a stable test particle's pericenter 



q = a{l — e) and apocenter Q = a{l -\- e) distances must 
obey 



q = a(l — e) > ai + ciRhi 
Q = a(l + e) < a2 - C2RH2, 



(6) 
(7) 



where the inner and outer planets are assumed to have 
circular orbits. We label the maximum stable eccentric- 
ity as etop (ffat top of the mountain) and consider the 
midpoint between the two planets (ai + a2)/2. We can 
rewrite Equations (|6| — ([7|) for particles on the edge of 
stability /instability as 



(ai +a2) 



(1 - etop) = tti + ciRhi (8) 



(ai +a2) 



(1 + etop) = a2- C2Rh2- (9) 



If we subtract the two equations from each other and 
solve for etop, we obtain 



Stop 



a2 — ai — CiRhi — C2RH2 



ai + a2 
CiRhi + C2RH2 
ai + a2 



a2 - ai 
ai + a2 



(10) 



We empirically determine ci and C2 by fitting Equa- 
tion (10) in a least-squares sense to values of etop mea- 
sured irom Figures [l]-|2] We find ci = 19.733 and C2 
= 4.1877. Comparison between computed values of etop 
(using Equation ([ip| )) and the measured values of etop 
(from Figures jl]-^) is shown in Figure [sj which illus- 
trates the maximum stable eccentricity etop as a function 
of planetary spacing for a range of planetary masses. 

The width of a stability mountain's base (where e = 0) 
can be related to the dynamical spacing criterion A be- 
tween two planets (see Equations (3| and Q). A system 
with two planets in a circular ana coplanar state satis- 
fies Hill stability (orbits do not cross) if A is greater than 
2>/3, or ~3.46 (Gladman 1993). The stability of systems 
with more than two planets are less well-characterized 
and are commonly determined using numerical calcula- 
tions. Estimates of the width (a2 — tti)stabie of each sta- 
bility mountain's base at e = (ignoring the effects of 
resonances as much as possible) are related to A (Fig- 
ure |4|. We do not find any stability regions for planetary 
systems with A < 10. 

We generalize the results shown in Figure [4] to a 
broader context. From this figure, we can determine a 
critical value of A that divides two-planet systems with 
stable versus no stable regions. This cross-over occurs 
in the range Acrit = 10 — 15. Accordingly, we suggest 
that two-planet systems similar to those explored in this 
paper cannot have extensive stability zones if their sep- 
arations have A less than 10. Similarly, we predict that 
stable regions can exist in systems with A greater than 
15. I n the February 2011 Kepler release (Borucki et al. 

)11 Lissauer et al.|20lTb| ), 95 out of a total of 115 two- 
planet systems have A > 15, or 82.6% of all two-planet 
systems in this sample can potentially harbor stability 
zones within the known planets. The results discussed 
here and illustrated in Figure [4l are consistent with previ- 
ous studies. ^ Chambers et al.^ (1996} numerically studied 
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Figure 3. Maximum eccentricity for a stable test particle 
orbit at a semi-major axis (ai + a2)/2 between two ex- 
isting planets. The unfilled squares represent estimates 
of etop with their uncertainties, and the filled circles and 
gray lines represent values of etop computed from Equa- 
tion ( [1Q| ) for various planetary systems (eight systems 
listed m Table [l] plus six additional systems for a larger 
sample). For comparison, the inverted triangle shows the 
planetary spacing between Jupiter and Saturn (note that 
we only considered two-planet systems with circular or- 
bits, and our results may not be applicable to systems 
with greater multiplicity of planets or non-circular or- 
bits). 
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Figure 4. Fraction of planetary separation (a2 —ai) with 
stable test particles at e = 0, as a function of planetary 
spacing criterion A. The unfilled squares represent the 
results from numerical simulations in this study (eight 
planetary systems from Table [l] plus nine additional sys- 
tems for a larger sample) with error bars representing 
measurement uncertainties. 

coplanar and circular configurations of 3, 5, 10, and 20 
planet systems, and found no stable systems with plane- 
tary spacing of A < 10. More recently. Smith & Lissauer] 
( |2QQ9 ) examined the packing density of systems with 3, 
5, and 9 Earth-mass planets in circular and coplanar or- 
bits with planets equally spaced in terms of A. They 
conducted long-term numerical integrations up to 10 bil- 
lion years, and demonstrated that 3-planet systems are 



stable when the spacing between neighboring planets is 
roughly A 7. Other previous results on spacing be- 
tween planets or protoplanets include the typical ~10 
Hill radii spacing between neighboring protoplanets, as 
seen in simulations of protoplanetary accretion from a 
swarm of planetesimals (Kokubo & I da 1998, 2000, 2QQ2| . 

5. SCOPE AND LIMITATIONS OF OUR RESULTS 

Given the large amount of possible parameter space 
that can be explored in stability studies, we summarize 
the limitations and scope of results stemming from this 
paper. We also discuss any other assumptions and con- 
siderations that may change our results. 

Sample. We solely investigated multi-planet systems 
announced by the Kepler t eam in the Februar y 2011 re- 
lease of candidate systems (Borucki et al.|2QlT]). No other 
planetary systems were considered. Therefore, our sam- 
ple has the same biases as any iTep/er detection, including 
the observational preference towards short-period plan- 
ets given the transit detection method. Our study is also 
limited to Kepler systems for which there are two known 
planets, and we do not investigate the dynamical spacing 
in systems with greater multiplicity of planets. 

Masses. The planetary masses are typically not known 
for these KOI systems. We estimated masses using 
Kepler-medisuied planetary sizes with a power law (Table 
[l] Lissauer et al.| ( |2Qllb| )) obtained from fitting to Earth 
and Saturn. However, the densities and true masses of 
Kepler planets can be different from these assumptions, 
which could change our results. 

Eccentricities. Eccentricity is another important dy- 
namical parameter that is not known for most Kepler 
multi-planet systems. We have assumed zero eccentrici- 
ties for the known planets in our numerical calculations, 
and this assumption is consistent with the expected tidal 
circularization of close-in planets. For the only confirmed 
planetary system in our sample, KOI 72 or Kepler-10, 
photometry and radial vel ocity data suggest th at Kepler- 
10b has zero eccentricity (iBatalha e t al.|2Q"TT ). Non-zero 
eccentricities of the planets in our sample, if present, 
would change the locations of stability regions of test 
particles. 

Inclinations. We assumed zero inclinations between 
the orbits of known planets in our sample as well as low 
inclinations up to '^5° for test particles. Consequently, 
our results can only be applied to systems that are rela- 
tively coplanar. The assumption of coplanarity or near- 
coplanarity is reasonable for multi-planet systems discov- 
ered by the transit technique at the heart of the Kepler 
mission, given that the inclination dispersion of these 
systems appears to have mean of < 10°. 

Integration time. We integrated test particles for a 
time span of 10^ years due to CPU time limitations, 
but more accurate modeling can be obtained by using 
a longer integration time period. There may be test par- 
ticles that are stable over 10^ years but not over longer 
timescales, although our simulations show that ~90% of 
particles unstable in 10*" years were unstable within the 
first 10^ years. 

6. CONCLUSION 

The "packed planetary systems" model advocates the 
idea that all planetary systems are formed to capacity. 
To test this hypothesis, we investigated the packing den- 
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sity of Kepler candidate two-planet systems from the first 
four and a half months of the mission. Through nu- 
merical calculations, we determined whether regions of 
stability exist between known planets with wide sepa- 
rations, i.e., in systems that seemed the most unpacked 
based on how well they satisfy Hill stability. Discovery of 
a stable region suggests that a low-mass body could be 
present in the gap, which would then bring the system 
to a more "packed" state. With time, such predictions 
will be shown to be correct or incorrect, allowing us to 
gauge the success of this model. 

We performed detailed numerical simulations of eight, 
two-planet Kepler systems, selected using an analyti- 
cal /3//3crit stability criterion. In addition to the known 
planets, we included 4000—8000 test particles per plane- 
tary system, allowing both circular and non-circular, and 
coplanar and non-coplanar orbits. These test particles 
are good proxies for low-mass bodies such as terrestrial 
planets as well as small bodies such as asteroids or dwarf 
planets. We integrated all bodies for 10^ years; we de- 
fined stable particles as those that remained stable dur- 
ing the length of the integration and unstable particles 
as particles that experienced a collision or ejection. 

Our results (Figures [l] to [2| indicated that all of the 
planetary systems investigated here (KOIs 433, 72, 555, 
1596, 904, 223, 1590, and 139) can pack additional, yet- 
undetected bodies in the identified stable locations. We 
also discussed relationships relating dynamical spacing 
between known planets and the extent of the inter-planet 
stability region. We derived an analytical relationship 
relating the largest possible eccentricity of a stable test 
particle to the semi-major axes and Hill radii of the two 
planets surrounding the particle. We also demonstrated 
that A, the separation between two planets in units of 
their mutual Hill radii, can be a reasonable predictor of 
whether or not stability regions can exist between plan- 
ets. The cut-off occurs at a critical A between 10 and 
15. We suggest that planets with separation A < 10 
are unlikely to host extensive stability regions. Based on 
this A = 10 — 15 cut-off, we suggest that about 95 out of 
a total of 115 two-planet systems in the February 2011 
Kepler sample may have sizeable stability regions. 

We thank the referee for useful comments. 
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